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We introduce a new numerical scheme for solving the initial value problem for quasiequilibrium 
binary neutron stars allowing for arbitrary spins. The coupled Einstein field equations and equations 
of relativistic hydrodynamics are solved in the Wilson-Mathews conformal thin sandwich formalism. 
We construct sequences of circular-orbit binaries of varying separation, keeping the rest mass and 
circulation constant along each sequence. Solutions are presented for configurations obeying an 
n = 1 polytropic equation of state and spinning parallel and antiparallel to the orbital angular 
momentum. We treat stars with moderate compaction ((m/R)^ = 0.14) and high compaction 
((m/R)oo = 0.19). For all but the highest circulation sequences, the spins of the neutron stars 
increase as the binary separation decreases. Our zero-circulation cases approximate irrotational 
sequences, for which the spin angular frequencies of the stars increases by 13% (11%) of the orbital 
frequency for (m/R) ao = 0.14 ((m/R)oo = 0.19) by the time the innermost circular orbit is reached. 
In addition to leaving an imprint on the inspiral gravitational waveform, this spin effect is measurable 
in the electromagnetic signal if one of the stars is a pulsar visible from Earth. 

PACS numbers: 04.30.Db, 04.25.Dm, 97.80.Fk 



I. INTRODUCTION 

Binary neutron stars are connected to two fascinating 
observational phenomena: 7-ray bursts and gravitational 
waves. Binary systems of compact objects (neutron stars 
and black holes) are strong emitters of gravitational radi- 
ation and thus prime candidates for the new generation 
of gravitational wave laser interferometers 0, II II 01 01 • 
Binary neutron stars may also be the engines for 7-ray 
bursts of short duration [6j. The numerical modeling of 
these systems becomes essential for the analysis of obser- 
vational data in both of these windows. 

Numerical simulations of binary neutron stars face two 
independent challenges: obtaining astrophysically realis- 
tic initial data and determining the numerical evolution 
of such data in time. In this article, we concentrate on 
the former (initial value) problem and extracting physical 
insight into the evolution of binaries from these solutions 



In thepast most researchers H El El El El El 
El El El ES El El El, have addressed the gravi- 
tational part of the initial value problem (IVP) for bi- 
nary neutron stars (i.e.; the solution of the Hamiltonian 
and momentum constraints for a quasi-equilibrium cir- 
cular orbit) via the Wilson-Mathews conformal "thin- 
sandwich" approach d H . This method is typically 
based on restricting the spatial metric tensor to a con- 
formally flat form and imposing a helical Killing vector 
to the spacetime to enforce the "circular orbit" condition 
[23|. The hydrodynamical part of the IVP however, has 
been addressed so far in three different ways: 

(1) Wilson, Mathews, and Marronetti jl used an 
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approximate set of relativistic Euler equations with post- 
Newtonian radiation reaction to relax the fluid motion 
to steady state. This method is difficult to implement 
numerically, since it requires a fully-functional, hydrody- 
namical evolution code. 

(2) Baumgarte et al. E3 solved the Bernoulli equation 
for a corotating binary. This corresponds to a static fluid 
in the reference frame rotating with the system, for which 
the Bernoulli equation applies. While straightforward to 
implement numerically, corotating systems do not rep- 
resent astrophysically realistic binaries, due to the high 
viscosity required to produce the tidal lock 0, |25j . 

(3) Bonazzola et al. EHHHJ, Teukolsky [H, and 
Shibata [2£j developed a formalism |3(j for systems with 
zero fluid vorticity (irrotational binaries). This formal- 
ism requires the solution of an extra elliptic equation for 
the matter (in addition to the Hamiltonian and momen- 
tum constraints). Adopting this method, Bonazzola et 
al. E3j Marronetti et al. |l(J, Uryu and Eriguchi [l7| . 
and Taniguchi and Gourgoulhon ; 21J obtained numerical 
solutions for the irrotational binaries. However, the extra 
elliptic equation requires boundary conditions that have 
to be specified on the stellar surface, something difficult 
to implement accurately in Cartesian coordinates. Dif- 
ferent investigators have resorted to different techniques 
to deal with this complication, like pseudo-spectral meth- 
ods 0, 0| , decomposition of the elliptic equation into 
homogeneous and inhomogeneous parts in Cartesian co- 
ordinates Il5l ITfl . and overlapping patches of spherical 
coordinates |l7lllc| . 

We present in this paper a new hydrodynamical recipe 
that has as its most important feature the capability of 
providing solutions to the IVP for binaries with stars with 
arbitrary spins. This scheme also requires the solution of 
an extra elliptic equation, but here the method reduces 
the difficulty attached to implementing stellar surface 
boundary conditions (particularly complex in Cartesian 
coordinates). This feature allowed us to obtain accurate 
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solutions using low resolution Cartesian grids. 

We construct sequences of quasi-equilibrium orbits 
with constant rest mass and constant relativistic circu- 
lation along the stellar equator. These sequences, which 
mimic inspiral evolutions outside the innermost stable 
circular orbit (ISCO), provide insight into the coupling 
of spin and orbital angular momenta in general relativis- 
tic binaries [3l| . 

In Sec. [H]we discuss our treatment of the gravitational 
fields, which follows the standard Wilson-Mathews con- 
formal thin-sandwich (CTS) approach. In Sec. II I II we 
introduce our new method that allows the free determina- 
tion of the stellar spins for configurations that satisfy the 
Euler and continuity equations. The full set of equations 
is summarized in Sec. Hvl Sections M and IvTl deal with 
our numerical implementation and code tests, including a 
comparison of our results with those of Baumgarte et al. 
[Tof for corotating binaries and Uryu and Eriguchi ^tJ for 
irrotational binaries. Finally, Sec. IVIII presents our re- 
sults for sequences of binaries in quasi-circular orbits for 
different values of the equatorial circulation, including a 
comparison with PPN results. 

II. SPACETIME EQUATIONS 

A. Initial value problem 

In this Sec. we present the equations for the Wilson- 
Mathews conformal "thin-sandwich" method [3, u| , fol- 
lowing the notation employed by Baumgarte et al. [Tof. 
We refer the reader to those papers for a more detailed 
explanation of the method and the justification of its ap- 
proximations. 

In the (3+1) spacetime formalism of Arnowitt-Deser- 
Misner (ADM) [3^, |3^| the line element may be written 

ds 2 = g^dx^dx" = - (a 2 - (3 t p l )dt 2 

+ 2j3idx l dt + -fi m dx l dx m , 

where g^ v is the four-metric, a the lapse function, (3i the 
shift vector, and 7^ is the spatial three-metric. Latin 
(Greek) indices go from 1 to 3 (0 to 3), and we set G = 
c = 1. The symbols ( ) ;(tt , Di, and V, will represent the 
covariant derivatives associated with the tensors <7 M „, 7^ , 
and Sij respectively. 

The Hamiltonian and momentum constraint equations 
can be written as 

R - K lm K lm + K 2 = 16np , 

D^RK-^K) = 8nf , (1) 

where R is the three-dimensional Ricci scalar, the 
extrinsic curvature, and K its trace. The source terms p 
and f are the total mass-energy density and the spatial 
components of the four-momentum density, and they are 
related to the matter stress-energy tensor by 

p = n^n u T^ v , 

f = -7> 7 T<"? , (2) 



where is the vector normal to a spatial hypersurface. 
Under maximal spatial slicing (K — 0), which we adopt, 
the third term and the second term of the left-hand side 
(LHS) of the Hamiltonian and momentum constraints 
vanish. We will also define for later use, a source term S 

S = j lm T lm . 

B. Reference frames and Killing vectors 

Binaries in quasi-equilibrium circular orbits are sys- 
tems that admit an approximate a timelike Killing vector 
of the form 

f = dt + n , 

in a coordinate frame tied to inertial observers at infinity 
at rest with respect to the system center of mass ("in- 
ertial frame"). Here f2 is the orbital angular velocity of 
the binary and <f> the azimuthal coordinate. In Cartesian 
coordinates, this vector has the components £ M = (l,£ l ) 
with = (Q, X rf (r the coordinate radius). No exact 
Killing vector of this type can exist in astrophysically re- 
alistic binaries, since the presence of gravitational radia- 
tion leads to an inspiral motion that breaks the symme- 
try. However the amount of energy dissipated per orbital 
period (outside the ISCO) is small enough to make this 
approximation valid in most scenarios [TlL 133. 13^1 13?3 | . 

In a rotating frame (i.e., a coordinate system that fol- 
lows the binary) the Killing vector adopts its simplest 
form £ = dt- Due to this simplification, we will work in 
such frame, where the stars are stationary j3]j • The rela- 
tion between the fields in the inertial and rotating frames 
is given in [3^ . 

C. Conformal thin-sandwich equations 

Following Q, the constraint equations |QJ are supple- 
mented by the time evolution equations for the spatial 
metric 

d tll3 = -2aKij + Difa + Djfa , (3) 

where Eq. J3J can be interpreted as the definition of 
the extrinsic curvature Kij . To simplify the problem the 
spatial metric 7^ is restricted to be conformally flat, 

la - . (4) 

where the conformal factor 4" is a positive scalar func- 
tion. For this conformally flat metric, the Ricci scalar R 
becomes 

R = -8*" 5 V 2 * , 

where V 2 is the flat space Laplacian. Using this equation 
and defining K l i = \p 10 AT 1 -' , the Hamiltonian constraint 
becomes 

V 2 * = --y- 7 K lm K lm -2tt* 5 p . (5) 
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To impose the condition of quasi-equilibrium, we set the 
time evolution of the trace-free part of the spatial metric 
7y to be zero. This restriction together with Eq. © and 
the definition of K l i lead to 



»!/'• / 2 , 

K 11 = — V l (3 J + V ] I3 1 - -S^Vtf 1 
2a \ 3 



(0) 



Now we can rewrite the momentum constraint equations 
as a set of elliptic equations for the shift vector 

1. 



V 2 /3 4 + -V^V^) 



2V; ln(a^- b )K u + IQiraWj 

These equations are split into two by decomposing (3 l 
according to 



p = - G l + (ft x r)* , 

with the corresponding equations for B and G l 
V 2 B = ViG 1 , 



V 2 G l = -2V;ln(a*- 6 )^ - IGna^f 



(8) 



Note that the term (fl x f) 1 that dominates the asymp- 
totic behavior of the shift vector in the rotating frame 
has been given explicitly, instead of including it in the 
fields G % . The reason is simply that this decomposition 
allows for better control of the boundary conditions for 
the Eqs. © (see Sec. HVIfor details). 

Finally, an elliptic equation for the product of the con- 
formal factor and the lapse function ('fa), is derived us- 
ing the time evolution equation for the extrinsic curva- 
ture Kij, by requiring the quasi-equilibrium condition 
d t K = 

V 2 (a*) = aV (^- s K lm K lm + 2irV 4 (p + 2Syj . (9) 



III. MATTER EQUATIONS 
A. Bernoulli equation 

We assume that the matter behaves like a perfect fluid 
with a stress-energy tensor 

= {p {l + e) + P) u»u v + Pg^ , 

with po the proper baryonic mass density, P the pres- 
sure, e the internal energy density, and m m the fluid four- 
velocity. As described in Sec. Ill Bl we will assume that 
the system is embedded in a four-dimensional manifold 
endowed with a Killing vector £ p . Following [3^|, we 
contract the Euler equations for a perfect fluid with the 
Killing vector £^ 

(p (l + e) + P)^«"« p .„ = -eP,» -eu»u v P, v (10) 



and derive the identity 



(ii) 



from the antisymmetry of £ p; „. Knowing that ^P,^ = 0, 
we combine Eqs. ffT77f) . (fTT|) 

(po(l + e) + P) d{u^)/dT = -u^dP/dr , 

which, together with the first law of thermodynamics 
d(p (l + e)) = (po(l + e) + P) dn /n , leads to 



dn _ rf(p (l + e) + P) 
n p {l + e)+P 



where n is the baryonic number density. Note that 
d (ln(tve) - In(no) + ln[p (l +e) + P]) = 0, 



or m more 



(7) compact form 



F = 



C n 



Po(l 



P 1 



(12) 



where F = and C is a constant. In general C is 

a constant only along flow lines and not all across the 
spatial slice. However, in two special binary systems, 
namely corotating and irrotational, C has been shown to 
be a spatial constant |2(J,|22|]. (In those two cases, in fact, 
C is a global spacetime constant since the fluid elements 
will carry the same value all along the flow lines.) Moti- 
vated by the fact that corotating and irrotational binaries 
describe extreme opposite physical systems, we will take 
as an approximation C to be constant across a spatial 
slice for our adopted velocity fields (see Sec. IIII C|l 40] . 

Now we will write the field F in a more tractable form. 
The Killing vector £ M can always be written as £ M = 
(!,£*)• Then F becomes 



F 



Uq + Ul£ 



(13) 



Substituting the contravariant components u l of the four- 
velocity 



Uq 

Uj 



9a ^ = goou 
3ip,y^ = g l0 u° - 



- goiu 1 
9iiu l , 



into Eq. @3J), and using Eq. (0J for the spatial metric, 
we get 

F = -a 2 u° + y% m (fl l + Z l )(u f3 m + u m ) . 

A more convenient form of this equation can be obtained 
by replacing the spatial components of the fluid four- 
velocity u 1 by the components of the fluid coordinate 
three- velocity v l = u 1 /u°: 

F = u° [-a 2 + y 4 6i m {/3 1 +£ l ){f3 m +v m )] . (14) 

The normalization condition for the four-velocity 



(15) 
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leads to the following expression for tt° 

u° = [a 2 - + v l )(f3 m + v m )] ~ 1/2 . (16) 

Inserting this expression into Eq. (|14fl we finally arrive 
at the formula 



F 



[a 2 - ^ 4 5i m {f3 l + v l ){f3 m + v m )} 



1/2 



(17) 



where F, and thus Eq. 1|12|) . is a function of the grav- 
itational and gauge fields a, and (3 % , Killing vector 
spatial components and the coordinate three- velocity 
of the fluid v l . 



B. Continuity equation 

In the previous Sec. we obtained the Bernoulli equa- 
tion l|12l) assuming the existence of a timelike Killing vec- 
tor and using the Euler equation for a perfect fluid. Now, 
we need to guarantee that the fluid will also satisfy bary- 
onic mass conservation 



(P0« M ) ;M = 



(18) 



In the rotating frame, Eq. (|18(l can be explicitly writ- 
ten for a spatially conformally flat spacetime as (see Ap- 
pendix ^ for the derivation) 



where 



Vi{ Po v l R ) + p v l R V l (ln(u°a^ 6 )) = , (19) 



v\ - (n x r) 1 , (20) 



and where R (I) labels a quantity in the rotating (inertial) 
reference frame. We can always decompose the three- 
velocity v R as the sum of a solenoidal and an irrotational 
part 



,J RS 



J RI i 



where 



and 



enm y' v^j = 



(21) 



(22) 



(23) 



and where eu m is the three-dimensional Levi-Civita ten- 
sor. The solenoidal component v RS will be specified to 
force the bulk of the fluid to move in a pre-determincd 
way (see Sec. IIII CI for more details). The irrotational 
component v RI on the other hand, will be determined by 
satisfying the continuity equation ltlD|) . In order to de- 
rive an equation for v RI , we will define a scalar function 
a such that 



The existence of this scalar is guaranteed by the condition 
(|23[l . Inserting Eq. 121|) into the continuity equation 
(|T9j> . we get 

Po^i(v RI ) +(v l R s + v l m ) 

(ViPo + Po V ; [ln( U a* 6 )]) - . 

We divide now by pa (which is non-zero inside the star) 
and collect the terms in the second parenthesis of the 
second term 

Vi(«k/) + («ks + «kr) V ; (ln( PoU °a* 6 )) = . 

Using the definition of Eq. (|24|l . we get an elliptic equa- 
tion for er 

V 2 cr = -{v l RS + V'cr) V/ (ln(p u°a* 6 )) . (25) 

This equation has to be solved in conjunction with an 
appropriate boundary condition. A typical choice is a 
condition that forces the fluid velocity at the stellar sur- 
face to be tan gent to this surface, like the ones employed 
in Newtonian |4l| and relativistic [2^| irrotational bina- 
ries. This condition can be written, for instance, as 

v 1 r VlPo\surf = (v RS + v/(J ) ^iPolsurf = . (26) 

A condition of this type becomes problematic in Carte- 
sian coordinates, since it involves imposing a Neumann- 
like condition on a near-spherical surface. Marronetti 
et al. 16] solved the problem by splitting the elliptic 
equation for the potential field into homogeneous and in- 
homogeneous parts, and letting the homogeneous field 
take care of the boundary condition. While formally cor- 
rect, this solution is very difficult to implement accu- 
rately and the results are strongly dependent on grid res- 
olutions 0. Because of boundary conditions like this, 
other groups resorted to spectral methods 0, |2(], 0] 
or superposing mappings of spherical polar coordinates 

003 

This problem can be minimized in the following way. 
Let us introduce the field A = crpo- Knowing that 

V/^oV'cr) = V 2 (crp ) - V ; (<tV%) , 

we can write the first term of Eq. (|19|) as 

Vi(A)t4) = Vi( Po v l RS )+Vi(p V l a) 



= Po^iv RS 



J RS 



v RI 



VV 



(24) 



+ V 2 (ap )- V ; ((tVVo) ■ 

Going back to Eq. (|19fl . using Eq. I|22|) . and rearranging 
the terms, we get an elliptic equation for the field A 

V 2 A - Vj(<7 VVo) - Pov l R V/ (\n(u°a^ 6 )) 

which, when writing the RHS in terms of A, becomes 

V 2 A = Vj(A V l ln(p ))-v l RS V^o+Viln(u°a* 6 ) 
(y l X-p v' RS - A Vjln(po)) • 
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To determine a boundary condition for this equation, we 
write the condition for a H26|) as a function of the new 
field A 

V'( — ) Vlp \surf = -Vr S V/p | Sur / . 

Po 

The RHS of this equation is zero if the stellar surface is 
spherical and the field v RS is specified by Eq. (|2~5)l given 
below. This will not be the case in reality, since tidal 
forces and spin will deform the surface away from spheric- 
ity. However, the spherical approximation was employed 
in previous work |15| . where it was shown to have little 
impact on the final results. Approximating the RHS of 
the above equation to zero, we can satisfy the resulting 
boundary condition by imposing 



A 



Surf 



= V'A 



Surf 



= o 



(27) 



in a vicinity of the stellar surface. For instance, this 
condition can be numerically imposed on the grid points 
where po falls below some threshold value (see Sec. HVjl . 
This approximate boundary condition is easier to imple- 
ment in Cartesian coordinates than the Neumann condi- 
tion JUl. 



C. Bulk fluid motion 

The bulk of the stellar fluid motion in the rotating 
frame resides in the field v RS . Since we are free to specify 
v RS in any way that satisfies Eq. 11221) , we will choose the 
form 



J RS 



(a — 1) d x (r — r ) 



(28) 



where a is a free parameter, Q is the angular three- 
velocity of the binary system, f is the coordinate position 
vector, and fo is the geometrical center of the star (deter- 
mined by the point of maximum baryonic density). This 
condition forces the bulk of the fluid to move like a rigid 
rotator with angular velocity parallel to the orbital an- 
gular velocity (l and magnitude (a — 1) fHn the rotating 
frame. The fluid three-velocity in the inertial frame v\ 
is related to its counterpart in the rotating frame v R by 
Eq. HI which, together with Eqs. and lf2~g|) 



gives 



v\ = (Q x r ) ,: + a [n x (r - Po)]* + VV 



(29) 



A star with this velocity profile will spin at a frequency 
(as seen by a distant observer in the inertial frame) Q s = 
a f2. Thus, we will refer to a as the spin parameter of 
the stars. If a = 1 then v RS is zero, and a — const is 
the exact solution to Eq. (|25|l , with boundary conditions 
(|26|l . Then v\ becomes (Q, x r) 1 , which is the condition 
that defines corotating binary systems. In appendix iBlwe 
show that for this particular case, the Bernoulli equation 
(l.')2t assumes the standard form for corotating binaries 



[see for instance Eq. (32) of Baumgarte et al. 0]]. Thus, 
setting a — 1 we recover the full set of CTS equations 
corresponding to corotating binaries. 

The spin parameter a allows us to control the stel- 
lar spin, while the scalar field a provides the correction 
necessary to satisfy baryonic mass conservation. Perfect 
fluid irrotational binaries have been obtained in the past 
by requiring the fluid vorticity tensor to be zero psl |29| 



(hu u );ft - {hu^.y = 



For these cases, the flow hu^ can be written as the gradi- 
ent of a potential field, for which a new elliptic equation 
is derived. The method we propose here does not enforce 
this condition on the fluid. Our method will rely on the 
concept of circulation to approximate solutions with zero 
vorticity. The Kelvin-Helmoltz theorem j4j| states that 
the rclativistic circulation 



C(c) = <j)^hu t _ l da f ' 



(30) 



is conserved for isentropic fluids following any closed 
path c, when evaluated on hypersurfaces of constant 
proper time. Integral conservation laws like the Kelvin- 
Helmholtz theorem can be defined more generally than 
on hypersurfaces of constant proper time. For instance, 
in a reference frame where exists a timelike Killing vec- 
tor whose coordinates are £^ = (1,0) , it can be shown 
than dC(c)/dt = [42j. This simply states the fact that 
in a frame where the fluid is stationary (like the frame 
corotating with the binary) the circulation is also time 
independent. Keeping the circulation constant from one 
orbit in the sequence to another is our working hypoth- 
esis. This is consistent with other assumptions of the 
quasi-equilibrium method, like assuming that the radial 
velocity of the stars remains zero all along the orbits of 
the sequence. In reality, the radial velocity like the circu- 
lation on t=const hypersurfaces will change slightly when 
the stars get closer to each other and full time numeri- 
cal evolutions will eventually estimate qualitatively this 
effect. In this paper we will evaluate the integral l|30|) 
along the stellar equator and this is what will be referred 
to as the circulation of the star. The two stars in the 
inspiraling binary will conserve the value of circulation 
Coo they had at large separation. The circulation of a 
single irrotational star (i.e., at rest in the inertial frame) 
is zero. Consequently, we will represent an irrotational 
sequence by the quasi-stationary CTS solutions for stars 
with C — 0. Note that conservation of circulation will 
also be true for sequences with C = Coo ^ 0. 

In the Newtonian limit, the value of a — will corre- 
spond to C = and in Appendix \C\ we show that in this 
case our scheme recovers the exact equations correspond- 
ing to Newtonian irrotational binaries. In Sec. IVIII we 
will show that, for the general relativistic case, our se- 
quences with C — agree very well with the irrotational 
binary sequences obtained using an exact treatment. Of 
course, by adopting the form of the velocity field accord- 
ing to Eq. H29f) we are making an approximation to the 
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exact irrotational velocity field. In the Newtonian limit, 
our velocity field requires a = to have zero vorticity. 

Finally, we note that Eq. Q29JI does not represent the 
only type of fluid motion that our method would allow. 
We could also specify a differential rotation velocity pro- 
file, by modifying the scheme presented in this paper to 
accommodate such circumstances. 



D. Equation of state 

The method described in the previous sections is in- 
dependent of the equation of state (EOS) used to model 
the stellar matter. In this paper and for the sake of sim- 
plicity, we will employ a polytropic EOS 



P = K p , 



(31) 



where k is the polytropic constant and T the adiabatic in- 
dex related to the polytropic index n by T = 1 + l/n. All 
the results presented here are for T = 2, which facilitates 
a comparison with previous work (see Sec. IVII|) . 

For polytropic EOS like Eq. the specific enthalpy 



h = exp 



dP 



Mi 



P 



becomes 



h = 1 + e H 



P 

A) 



Note that this expression, combined with Eq. (|12H pro- 
vides a very compact form for the Bernoulli equation 

F=j, (32) 

where C is a new constant 01 ■ 

It will be convenient to define a variable q = po/P and 
write every quantity as function of q [T(i| 



Po 
h 
e 
P 



K- n q n , 

q (n + 1) + 1 

n q , 

K - n q n+1 . 



The matter source terms can then be written as 
p = po (l + «(l + n)) (au°) 2 -P , 
f = po (1 + q(l + n)) u° 2 a (v i R + p l ) , 
S = 3P + po (1 + q{l + n)) [(cm ) 2 - 1] . 

Using these relations and the normalization condition 
for the fluid four-velocity i|15|) we can derive the expres- 
sions needed for the RHS of the elliptic equations JSJ, 
©, and @ 

p = K - n q n [(l + q(l + n)) (au°) 2 -q] , 

f = K - n q n (l + q (l+n))au o2 (v* R + P) , 
p + 2S = K - n q n 

[(l + g(l + n)) (3 (cm ) 2 - 2) + 5g] . 



TABLE I: Boundary conditions for the outer boundaries 
(r — > oo) and on the coordinate planes in Cartesian coordi- 
nates. The stars are aligned along the y axis with a coordinate 
separation of d and the equatorial plane is z = 0. 



r —* oo 


x = 


y = o 


z = 


* - 1 ~ f(r) 


even 


even 


even 


(a*) - 1 ~ f(r) 


even 


even 


even 


/-ix y 


even 


odd 


even 


T 

— 3 


odd 


even 


even 


z xyz 


odd 


odd 


odd 


B ~ ^ 


odd 


odd 


even 



^y 2 



IV. FINAL SET OF EQUATIONS 



In order to simplify the numerical problem, we find 
useful to nondimensionalize all the equations by factor- 
ing out the polytropic constant k. Following ^(ji we 
use dimensionless coordinates x l — n~ n / 2 x l , derivative 
operators V 1 = «; n ' 2 V*, masses M — n~ n / 2 M, angular 
velocities fl = K n / 2 £l, angular momenta J — n~ n J, etc. 
This is equivalent to assuming a value of K — 1 and omit- 
ting the bars on top of the new variables and operators. 
In the remainder of the paper, all the variables and quan- 
tities will be given in these dimensionless units. Results 
for arbitrary values of k can be determined by applying 
the above scaling relations. 



In the previous sections we derived a set of CTS equa- 
tions that define our binary system. The fields that will 
play the role of independent variables are the conformal 
factor 'J, the product of the conformal factor and the 
lapse function (a^>), the vector field G % and the scalar 
B, from which the shift vector j3 l is derived, the flow 
potential A and the matter density variable q. The cor- 
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responding equations are 
1 



V 2 (a#) 



V 2 G l = 



2tt4'V 



[(l + ?(l + n)) (c^ ) 2 -?] , 
8 

[(l + g(l+n)) (3(cm ) 2 - 2) + 5 q] 
-2Vi\n{a^- & )K a - 16tt*V 
(4 + ^) (l + «(l + n)) (mi ) 2 , 



V 2 B = V Z G' , 

V 2 A = nVi{\W l hx{q))-^ RS Wiq n 
+ (V'A - q n v l RS - A n V 1 ln(q)) 
Vzln(u°a* 6 ) , 



and the Bernoulli equation 



n+ 1 V F 



g-1 



(33) 



(34) 
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FIG. 1: Hamiltonian constraint residual across one of the 
stars. The residual is evaluated along the line (0.22, y, 0.22). 
This line is chosen off-centered to avoid the cancellation of 
some of the constraints due to the symmetries of the prob- 
lem. The residuals have been plotted for three different grid 
resolutions: 32 3 (dotted line), 64 3 (dashed line), and 128 3 
(solid line). 



These equations are completed with the definition of 
©, the formulas that define the relation between the 
shift vector f3\ G\ and B Q, the field F <[T?jl. and the 
velocity field v RS P5)l. 

We work with binary systems composed of identical 
stars aligned along the y axis with z = as the orbital 
plane. The symmetries of these systems allow us to solve 
the equations in one octant of the numerical grid. Equa- 
tions l|33|) form a set of seven nonlinear, coupled elliptic 
equations that is solved numerically as described in Sec- 
tion The boundary conditions for the first six equa- 
tions are specified at large distances (i.e., at the outer 
boundary of the numerical grid). They follow a Robin- 
like algorithm based on the fall-off behavior of the lowest 
order non-vanishing multipole moments for the fields G l 
and B, and the first and second non- vanishing multipole 
moments for the fields VP and (a\&). The quadrupole mo- 
ments for VP and (a^) are required in order to get an 
accurate measure of the gravitational (ADM) mass, due 
to the close proximity of the boundaries. The details of 
the fall-off dependence for each field, as well as the re- 
flection symmetries at each plane, are given in Table [I] 
The elliptic equation for the field A requires A = V'A = 
at the surface of the star. This condition is imposed by 
forcing A = whenever the baryonic density falls below 
some threshold value, which is usually set at 5% of the 
maximum stellar density, to avoid problems related to 
the steep gradients of the density right at the stellar sur- 
face. This threshold can be dropped to smaller values by 
increasing the grid resolution. 

In Appcndix[D]we derive the formulas used to calculate 
the rest mass Mo, the ADM mass M, the angular mo- 
mentum J, and the circulation C of each quasi-stationary 
orbit. 



V. NUMERICAL IMPLEMENTATION 

A. Elliptic solver 

The core of the numerical code used to solve the set of 
equations (|33|l is an elliptic solver. This solver is based 
on a variation of a finite difference multigrid algorithm 
like the one used for previous work |13l ll(il| . The solver 
makes use of the octant symmetry of the problem and 
works in the Sec. of the physical space with positive 
values for the coordinates x, y, and z. Reflection sym- 
metries and boundaries conditions are described in Sec. 
IIVI For the results presented in this article, three differ- 
ent grid sizes were employed: small (32 3 ), medium (64 3 ), 
and large (128 3 grid points). Larger grid sizes might be 
needed for the generation of initial data sets for evolu- 
tionary codes and will be discussed in future work |43| . 
A more detailed description of the iterative algorithm at 
the core of the solver is provided in Appendix [EJ Com- 
puting the solution of Eqs. (|33|l for a single separation 
distance on the large grid takes typically around 50 CPU 
hours on the IBM Regatta p690 located at the National 
Computational Science Alliance (NCSA). 

VI. CODE TESTS 

We tested the code solving the set of equations (I33|l and 
the Bernoulli equation (|34|l for a binary with identical 
stars with individual rest mass mo = 0.1469, total mass- 
ener gy m = 0.1368, and compaction ratio (m/R)^ 
0.1 I |4J]. The stars were set along the y axis, at a coor- 
dinate separation between centers of d — 3.50, where the 
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FIG. 2: Momentum constraint residuals across one of the 
stars. The residuals are evaluated along the same line as in 
Fig. Line captions are identical to those of Fig. 



TABLE II: Comparison between the correct numerical solu- 
tion for a binary with stellar rest mass mo = 0.1469, coor- 
dinate separation d — 4.12, and circulation C = 0, and the 
approximate solution obtained when the flow potential A is 
neglected. As expected, the influence of the field A in the 
final results is small. 



Parameter Correct Result A = Result Diff. (%) 



Q, mo 

M 

J 

Max. po 
a 



0.00858 
0.27174 
0.08494 
0.12375 
0.58 



0.00843 
0.27160 
0.08268 
0.12381 
0.60 



1.7 

< 0.1 
2.7 

< 0.1 
3.3 



0.0015 



0.001 






(32) 3 x1 

(64) 3 x4 

(128) 3 x16 


0.0005 


n nnnc 















FIG. 3: Continuity equation residual across one of the stars. 
The residuals are evaluated along the same line as in Fig. 
Line captions are identical to those of Fig. Q 



stellar center is defined as the point of maximum bary- 
onic mass density po. The system rotates counterclock- 
wise in the orbital plane z — and the spin parameter a 
is set to be zero. This test problem was solved for three 
different grid sizes with 32 3 , 64 3 , and 128 3 grid points 
respectively. These grid sizes correspond to resolutions 
of about 15, 30, and 60 grid zones across the stellar di- 
ameter. In order to check the second order convergence 
of the code, we plotted the residuals of the Hamiltonian 
constraint (Fig. 0, the components of the momentum 
constraint (Fig. |2J), and the continuity equation (Fig. |3J), 
all evaluated in the stellar interior. The residuals are de- 
fined as the absolute value of the RHS minus the LHS of 
Eqs. and l|18fl respectively. The plots show their val- 
ues along an off-centered line defined by the coordinates 
(0.22, y, 0.22) that crosses the upper star from side to 
side. This off-centering was done to prevent cancellations 



in the constraints due to the symmetry of the system. 

Figure 0] shows the contour plot of the baryonic den- 
sity and the fluid three- velocity v R in the rotating frame, 
while Fig. 03 shows the irrotational component of the ve- 
locity v l RI (here Mq = 2mo). As we anticipated in Sec. 
IIII CI most of the fluid motion follows the solenoidal field 
v RS , determined by Eq. 1)28(1 . By setting v RS , we fixed 
the dominant part of the fluid velocity, while letting the 
field v l RI play the role of small corrections that guaran- 
tee satisfaction of the continuity equation l|18f) . Note that 
Fig. also highlights another convenient feature of our 
method: the fact that the field A, which is the most dif- 
ficult numerically, contributes minimally to the final re- 
sults for the cases explored here. While the velocity field 
v RI is at its largest near the surface, the matter current 
Po v ri is minimal due to the very steep decay of the rest 
mass density with increasing radial distance (see contour 
lines in Fig. [SJ). To verify that only small fractions of 
the total mass and angular momentum of the system are 
associated with the field v RI , we compare the main pa- 
rameters of two identical orbits, one of which was solved 
neglecting completely the flow potential A (i.e., A = ev- 
erywhere) . The difference between the correct numerical 
results and the ones without A is of the order of a few 
percent at most. This result justifies our approximate 
boundary condition for the scalar potential A l|27|) (see 
Table HTJ) . 

Numerically, our solutions obtained in Cartesian coor- 
dinates are as accurate as the more sophisticated (but 
more complicated) alternative numerical methods that 
employ patches of spherical coordinates These so- 

lutions still do not reach the very high accuracy of the 
pseudo-spectral methods [12, W$ , but Cartesian coordi- 
nates methods usually behave better in the presence of 
cusps at the stellar surfaces (when the stars are very close 
to each other). Cartesian coordinates are also more con- 
venient when providing initial data for evolutionary codes 
that work on those coordinates, since interpolating the 





X/M n 





X/M 



FIG. 4: Contour plot of the full fluid (coordinate) three- 
velocity v R and the rest-mass density po in the rotating frame. 
The contour lines are drawn for p = 10" (0 ' 2 3+0 A) po' ax , 
where po Iax denotes the maximum value of the rest-mass den- 
sity po (here it is 0.1241), for j — 0,1..., 7. Vectors indicate 
the local velocity field and the scale is as shown in the top 
right corner. The coordinates have been normalized by the 
total rest mass (Mo = 0.2938). 



fields from one coordinate system to another usually in- 
troduces noise into the data set ^3 ■ Figure shows 
the rest-mass contour plot for the closest orbit of the se- 
quence (™/-ff)oo = 0.14 and C = |45(. No equilibrium 
solutions were found inside this orbit. 

Friedman et al. [46( derived a virial relation for 
quasi-equilibrium spacetimes like the ones described here, 
showing that the virial theorem is satisfied if the grav- 
itational mass of the system is identical to the Komar 
mass. In the solutions presented in the next section, the 
relative difference between these masses is of about 2%. 
This difference is due to the proximity of the grid outer 
boundaries, since both masses are only defined asymp- 
totically. Test runs performed with various grid sizes 
show convergence of the relative difference to zero. We 
also find that the quasi-equilibrium relation Q — dM/dJ 
along the equilibrium sequences is satisfied to within 8%. 



VII. RESULTS 

All the binary sequences described in this paper where 
modeled using a polytropic EOS with index n=l (r = 2). 
For this particular EOS, the critical rest mass (grav- 
itational mass) of a star in isolation is mo = 0.180 
(m = 0.164) with a compaction ratio of (m/R)^ = 0.216 
[47| . All the binaries are composed of identical stars and 
the orbital plane is z = 0. 

We study models that have two different compaction 



FIG. 5: Contour plot of the irrotational component v RI of 
the fluid coordinate three-velocity and the rest-mass density 
in the rotating reference frame. The plot is labeled as in Fig. 
□ 
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FIG. 6: Contour plot of the rest-mass density po for the 
closest orbit of the (m/R)^ = 0.14 and C — sequence. The 
plot is labeled as in Fig. 0] 



ratios in isolation: a moderate value (m/R)^ = 0.14 
and a high value (m/R) x — 0.19. These compaction 
ratios correspond to individual stars with rest masses 
mo = 0.1469 and mo = 0.1767, respectively. When 
these stars are spun up to the mass-shedding frequency, 
their equatorial circulation assumes the critical values 
C C rit = 1-96 and C cr i t = 2.02, respectively, values which 
will be used to normalize the circulations of the binary 
sequences in the next section. 



10 




fl m Coordinate Separation 



FIG. 7: Angular momentum vs orbital frequency. We com- 
pare a sequence of corotating orbits for stars with compaction 
ratio at infinity (m/R)^ = 0.15 from Baumgarte et al. |To|| 
(dashed line) with our sequence with a = 1.0 (filled circles). 
For the irrotational case, we compare a sequence of orbits for 
stars with compaction ratio at infinity (m/R)^ — 0.14 from 
Uryvl and Eriguchi flTI) (solid line) with our sequence with 
null circulation (filled squares). 
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FIG. 8: Binding energy comparison between our results and 
previous work. The binding energy is defined as BE = (M — 
Moo) /Mo, where M is the gravitational (ADM) mass. The 
notation is identical to that of Fig. Q 



A. Comparison with previous work 

We tested our code by reproducing a corotating se- 
quence (a = 1) identical to one provided by Baumgarte et 
al. ^3 for stars with compaction ratio (m/R)oo = 0.15. 
We also compared one of our C — sequences with the 
irrotational counterpart given by Uryu and Eriguchi 0] 
for stars with compaction ratio (m/R)oo — 0.14. Figures 
and |H1 show the total angular momentum and binding 
energy versus the orbital angular velocity. The agree- 
ment between our a — 1 sequence and the corotating 



FIG. 9: Ratio of spin vs orbital angular velocities as a func- 
tion of the coordinate separation. The filled circles (open 
squares) correspond to the (m/R)^ = 0.14 (0.19) and C = 
sequence. The solid lines show the PPN prediction (Appendix 
EJ. 



sequence from 

.10] 

further validates our numerical code. 
Most significantly is the close agreement between our 
C = sequence and the irrotational sequence from |17| . 
given that the two sequences were obtained using dif- 
ferent (albeit closely related) formalisms and numerical 
approaches. The overlapping system of spherical coor- 
dinates used in 0] is better suited than the Cartesian 
grids used here for dealing with the boundary conditions 
of the elliptic equation for the fluid velocity potential. 
However, our splitting of the velocity field into solenoidal 
and irrotational parts compensates for this handicap to 
the extent that we obtain results of similar quality to 
those of 0] with modest grid resolutions (from 20 to 30 
grid points across the star). 

An interesting effect is revealed on the plot of the spin 
parameter a = f2 s /f2 versus coordinate separation be- 
tween stellar centers. Figure [5] shows the spin parame- 
ter for the C = sequences for both compaction ratios 
(m/R)oo = 0.14 and (m/R)^ = 0.19. We see that the 
stellar spins, which point along the direction of the or- 
bital angular momentum, increase as the stars get nearer. 
For comparison we also show the corresponding post- 
Newtonian prediction. As it can be seen, the agreement 
is close all the way up to the closest orbits (see Appendix 
|0 for more details and post-Newtonian derivation). 

B. Constant circulation sequences 

We constructed five sequences for the 
(m/i?)oo = 0.14 binary, with circulations 
C/C crU = (-0.510,-0.255, 0.0, 0.255, 0.510), and 
three for the (m/R)oo — 0.19 system with circulations 
C/C cnt = (-0.124, 0.0, 0.247). " 

Figures 1101 and 1111 show the angular momentum and 
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FIG. 10: Angular momentum (upper) and binding energy 
(lower) vs orbital angular frequency for the sequences with 
(m/i?) oo = 0.14. The lines are polynomial fits to the data and 
the labels denote the circulation ratio C/C cr it for the curves 
from top to bottom (upper plot). 



FIG. 11: Angular momentum (upper) and binding energy 
(lower) vs orbital angular frequency for the sequences with 
(m/R)oc = 0.19. The labels denote the circulation ratio 
C/Ccrit for the curves from top to bottom (upper plot). 



binding energy as a function of the orbital frequency for 
the cases (m/R)oc = 0.14 and (m/R)oo = 0.19, respec- 
tively. As expected, the angular momentum increases 
with the circulation C, due to the addition (or subtrac- 
tion, when C < 0) of spin angular momentum to the 
system. In the same way, the binding energy increases 
with \C\ due to the extra rotational kinetic energy. Note 
that none of these curves show turning points prior to ter- 
mination, in contrast to corotating sequences |10| . With 
the exception of the (m/R)oo = 0.14 and C = 0.0 se- 
quence (Fig. 0), we did not continue the orbits all the 
way up to the point of contact for two reasons: first, the 
code requires very high resolution to resolve the cusp at 
the contact point, making the calculation quite resource- 
intensive. Second, the closest orbits treated here appear 
to be inside the ISCO, thus rendering any analysis of the 
physical aspects of such orbits meaningless 48]. This 
claim is based on the dynamical determination of the 
ISCO that will be reported in a forthcoming paper [43| . 
in which the initial data sets presented here are evolved 
for a couple of orbits. Also, due to the high resolution 
and large CPU time requirements, only three sequences 
are presented for the case (m/R)^ = 0.19. Future work 
contemplates a more extensive cov erag e of the range of 
compaction ratios and circulations |49j. 

Figure ^] shows the spin vs orbital angular velocity 
for both sets of sequences, while Fig. ^| compares only 
the C = sequences for both compaction ratio cases. 
The neutron stars spin up during the inspiral, starting 
from zero spin at infinity to spin frequencies that are 
13% (11%) of the orbital angular frequencies for the se- 
quence with (m/R)oo = 0.14 (0.19). This spin-up effect 
is observed in all the sequences with the exception of the 
(m/i?)oo = 0.14 and C = 0.510 case (the largest circu- 
lation value studied in this paper), where the spin rate 
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FIG. 12: Spin vs orbital angular velocities for the sequences 
with (m/R)^ = 0.14 (upper) and (m/R) ao = 0.19 (lower). 
The line captions are those of Fig. ^] (upper) and Fig. fTTI 
(lower). 



decreases. This spin-up effect is imprinted in the inspiral 
gravitational radiation waveform and potentially observ- 
able in the electromagnetic spectrum if one of the neutron 
stars is also a pulsar (radio or X-ray). 

Figure 1141 shows the ratio of spin angular velocity vs 
orbital angular velocity as a function of the coordinate 
separation. 



VIII. CONCLUSIONS 

We introduce a new method for constructing quasi- 
equilibrium sequences of binary neutron stars in circu- 
lar orbit, allowing for the free specification of the rota- 
tional component of the fluid velocity. This allows us 
to choose the spin of the stars, thus extending the tra- 
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FIG. 13: Spin vs orbital angular velocities for the C = 
circulation sequence for compaction ratios (m/R)oa = 0.14 
(solid line) and (m/R)^ — 0.19 (dashed line). 
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the quasi-equilibrium binary formalism). In particular, 
we find that the nonlinear gravitational fields induce a 
spin-up in the stars for all the cases, except the sequence 
with the largest value of C. The variation of the spin of 
the star throughout the inspiral is a potentially observ- 
able electromagnetic effect when one of the neutron star 
is a pulsar. 

We plan a more extensive study of the parameter space 
(Mq/C) for different EOSs and a detailed study of the be- 
havior of the central density vs separation throughout the 
inspiral |49j| . The quasi-equilibrium orbits presented here 
are particularly well suited as initial data sets for evolu- 
tionary codes. In a forthcoming paper, we will study the 
evolution of such a set for binaries at different separations 
tracking them for the last couple of orbits to determine 
dynamically the location of the ISCO [43j . 
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APPENDIX A: CONTINUITY EQUATION IN 
CONFORM ALLY FLAT SPACE 



FIG. 14: Ratio of spin vs orbital angular velocities for the 
sequences with (m/R) ao = 0.14 (upper) and (m/R)oa = 0.19 
(lower). The line captions are those of Fig. 1101 (upper) and 
Fig. EI] (lower). 

ditional options of corotating and irrotational binaries. 
The sequences of circular orbits are characterized by two 
free parameters: the baryonic mass of each star, as in the 
standard quasi-equilibrium approach, and the relativistic 
equatorial circulation C. The conservation of the latter 
during inspiral along hypersurfaces of constant proper 
time is guaranteed by the Kelvin-Helmoltz theorem |42fl . 
Additionally, the numerical formalism is designed to pro- 
duce high quality initial data sets using Cartesian grids 
with modest resolution for future evolution studies. 

We constructed sequences for different values of the 
circulation for binaries with stars of moderate and high 
compaction ratio. Each circulation value corresponds to 
a particular value of spin (perpendicular to the orbital 
plane) when the stars are at infinite separation. These 
sequences show the spin-orbital angular momentum cou- 
pling that accompanies the inspiral (in the framework of 



The continuity equation l|18|) expresses the conserva- 
tion of baryonic mass in the system. In this appendix we 
will derive the expression for such equation in the rotat- 
ing frame when the spatial metric tensor is conformally 
flat. 

The continuity equation can be split as 

( P 0<);m = V>oO + Po T^u R = . (Al) 

Since the system is stationary in the rotating frame, any 
partial time derivatives vanish, reducing the first term of 
the RHS of Eq. (MJ to 

V m (po<) = VoGooO + Vi{p u l R ) = Vi{p u l R ) . (A2) 

We concentrate now on the term po r^„u^. Using the 
definition of the affine connection, we see that the con- 
traction r^ vanishes in the rotating frame 

= ^ CT (v^ CTO -v a5M) ) = o. 
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The term po T^ u u R gets reduced to the sum over the 
spatial components of the four-velocity u R 

Po = po r^u l R . (A3) 

In order to evaluate the sum over the affine connection 
components we will use the conformally flat spatial 
metric introduced by Eq. (J2J , together with the relations 
between the lapse, shift vector, and spatial metric and 
the components of the four dimensional metric tensor 
9fiv H3- After some algebraic steps, we get 

rg, = - V ia - -L [* 4 (/3 m V ; /3 m + /3 m V m/ 9') 
a 2a z 

+ 4$> 3 /3 l j3 m V m *] 



6_ . 1 
4* 3 /3 z /3 m V m *] 



where the raising (lowering) of indices is done using 8 lJ 
(5ij) which we have omitted for simplicity. With this 
result, we write Eq. (| A3|) as 

po t »v u r = po rj>k 



= p u R -V t a+ -V/* 
= Po u' R V ; ln(a* 6 ) . (A4) 
Replacing Eqs. (|A2|) and ((A4(l in Eq. (|A1|) . we get 

Vi(/5 ?4) + i oo^Vi[ m («* 6 )] = . 
which we divide by u° 

^i(pou l R ) + po^V/Iln^ 6 )] = . (A5) 
The first term can be decomposed 

= ^(.o^+Po^V.dn^)), 

which, when combined with the second term, leads to a 
new expression for the continuity equation (|A5|) 



For convenience, we will use the fluid three- velocity v l R = 
u R /vP as our variable, thus arriving at the expression 

V;(Po*4) +Po^Vi(ln(ii°a* 6 )) = . 



APPENDIX B: BERNOULLI EQUATION FOR 
COROTATING BINARIES 



Setting a = 1 in Eq. (21 gives the three-dimensional 
coordinate velocity in the inertial v\ = (fi X r) 1 , which is 
identical to the spatial components of the Killing vector 
£j in the same frame. We can write the Killing vector as 

^ = (Mi) = (i,4/u°) . 

As expected, the Killing vector is proportional to the 
fluid four- velocity: £j = /u°. The field F becomes 

F = uj^ = -^ , 

where we have used the normalization condition 1)15(1. 
Using this expression together with Eq. l(52*(l . we get 



— = const , 
h 

which is the familiar expression of the Bernoulli equa- 
tion for corotating binaries with polytropic EOS (see for 
instance Eq. (32) of 



APPENDIX C: NEWTONIAN LIMITS 

We derive in this appendix the Newtonian limit of the 
hydrodynamical equations employed in our method. The 
gravitational part is simply given by Poisson's equation 
for the Newtonian potential tp^. We start by keeping 
only the terms to first order in the Newtonian potential 
4>n, the Newtonian enthalpy h^, the rest mass density 
POi the square of the fluid velocity in the inertial frame vj, 
and the contraction 8i m (Cl x rfvf 1 - The simplest limits 
are 



h 

a 



1 + h N 
($1 x f)* 



These limits in combination with the equation for the 
time component of the fluid four- velocity u° give 



1 - (t> N + vj/2 , 



(CI) 



which together with — (SI x r) 1 , yield the limit of the 
field F jn|l 

F —> -[1 + 4>N + vj/2 - 6 lm {n x r) l vf] . 

Now the limit of the Bernoulli integral 1(32(1 is simply 

1 + <j) N + vj/2 - s lm (n x rfvf +h N = C . 



Making C = C — 1 a new constant, we get the Newtonian 
Bernoulli equation 



to + vj/2 - 5 lm {il x r) l v? + h 



N 



c 



(C2) 
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The Newtonian limit of the elliptic equation for the 
potential a is easy to obtain by realizing that the limit of 
the argument of the logarithm in the RHS of Eq. (|25|) is 
just po. Then, the Newtonian limit of Eq. (|25[) becomes 

VV = -(^5 + VM^- (C3) 

In similar manner, we get the limit for the elliptic equa- 
tion for the flow potential A 

V 2 A = V,(A V* ln(p )) - v l RS V/p . (C4) 

Equations l(U2)l and (|U3|l [or equivalently Eq. (EH)], 
together with Poisson's equation for the Newtonian po- 
tential cf) _/y form the complete self-consistent set of equa- 
tions necessary for applying our method in Newtonian 
gravitation, appropriate for scenarios like binary white 
dwarfs, etc. Note that in the case C = and a = 0, 
Eqs. I|C2|) and (|C3ll reduce to the equations derived by 
Bonazzola et al. [41| for Newtonian irrotational binaries. 

APPENDIX D: CALCULATION OF MASSES, 
ANGULAR MOMENTUM, AND CIRCULATION 

Every quasi-stable circular orbit in each sequence is 
characterized by a set of global quantities: the rest mass 
Mo, the gravitational or ADM mass M, the total angular 
momentum J, and the circulation along the stellar equa- 
tor C. In this appendix we present the formulas used to 
evaluate those quantities using the "k = 1" unit system 
described in Sec. IIVI 

The rest mass Mq of the system is [38| 

M = [ t> 6 au°q n d 3 x , 

where we have used po — q n . The total gravitational 
mass of the system can be read off of the monopole term 
of the conformal factor 'J, given that if? = 1 + M/2r + 

M = — - I \7 i ^dS l = -— [ V 2 fd 3 x . 

2^ Too 27T 

This last volume integral is evaluated replacing V 2 ^ by 
the RHS of the first (Hamiltonian) equation of the set 
(El- 

The orbital angular momentum of the binary is (by 
construction) positive along the z direction. Since we 
only consider stars with spins along the z axis, the total 
ADM angular momentum of the system remains confined 
to this direction. Following |5l| we write 

j = Jz = e -^tl x >K lk dSi 

= [ x^V l K lk d 3 x . 

8tt Joe 



Using the identity Xl\K lh = <f? 10 DiK and the momen- 
tum constraint Q we get 

J = I * 10 (x j y - y f) d 3 x 

J OO 

= f * 10 q n [1 + (1 + n)q] u° 2 a 

J OO 

[x (v y R + P y )-y (v x R + n] d 3 x, 

where we have used some of the expressions derived in 
Sec. liTTDl 

The relativistic circulation C, defined by Eq. (|3(J|) 
presents the numerical complication of requiring a line 
integral along the equator which, being a quasi-circular 
path, makes the result sharply dependent on the resolu- 
tion of our Cartesian grid. To minimize this effect, we 
transform the line integral into a surface integral using 
Stokes theorem 

[ dQ= [ 9 . (Dl) 
J o Jan 

In the case of Eq. (|3U[I . we have 
6 = hu^dx^ , 

which gives 

de = d(h u^) A dx^ + h A d{dx tl ) 
= (h u^), v dx v A dx^ . 

Using Eq. IjDljl . we can write the circulation as 

C(c) = <j> hu^da^ = J [(h u y ) x - (h u x ) y ]dS , 

where 5* is the surface confined by the equator on the 
z = plane. 



APPENDIX E: ELLIPTIC SOLVER ITERATIVE 
ALGORITHM 

The set of elliptic equations (|33fl and the Bernoulli 
equation (|34fl are solved iteratively, starting with an ini- 
tial guess for the fields. The iteration process is illus- 
trated by the following pseudo-code: 

Set Binary parameters : mo (= Mo/2), d, Coo 
Set Initial Iteration Parameters : £1 , (5£l)o 
Set Fields Initial Guess 

Do iter 
Do index = 1,7 

olve(V t lter - p lter ) 
End Do index 
Do j = -1, 1 

= Qiter-l + j {Styiter-1 

Compute q~j using (lj 
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\Qj — 1iter-l\ 1 2 



min((Ag)_ 1 ,(Ag) ,(Ag) 1 ) 



rpindex I 
r iter-l\\2 



Calculate (Aq)j 
End Do j 

J = j, such that (Aq)j 
Set q it er = qj 

Set Qiter = ^1,7 

If (J = 0) then 

(SQ)iter = (^)iter-l/ 2 

else 

(£f2)iter = (<551)iter-l 

End If 

Find C such that the rest mass is mo 
Find a such that the circulation is Coo 
Compute (AF) mdex 
If (AF) mdex < e) for (1 < index < 7) then STOP 
End Do iter 



The following is a brief description of the main features 
of each code section: 

Set Binary Parameters: The three free parameters of 
an identical star binary system are set. For our 
code we have chosen the individual star rest mass 
mo, the coordinate separation distance d, and the 
equatorial circulation at infinite separation Coo ■ 

Set Initial Iteration Parameters: We start the itera- 
tion providing an initial value to the binary angular 
velocity Ho and its error bracket (<5f2)o- This last 
value will be reduced during the iteration process, 
until it reaches the desired accuracy. 

Set Fields Initial Guess: The gravitational fields and 
the rest mass density q are initially set by mapping 
the solution of the TOV equations onto the octant 
grid. 

Loop iter: Main relaxation loop of the code. 

Loop index: This loops goes across the list of seven el- 
liptic equations. The source term and boundary 
conditions of equation index are computed using 
the values of fields 1 to (index — 1) at iteration step 
iter and the values of fields index to 7 at iteration 
step (iter — 1). 

Loop j: This loops determines the next value of the 
angular velocity fi. The rest-mass density field 
qj is computed using the Bernoulli equation 
for three different values of the angular velocity 
flj, namely (% ter _i) -(6D.)^ ter _^), fi(; teT -i), and 
(0( iter _i) + (50)( iter _i)). The value of fi; ter is set 
to the value of VLj that minimizes the Li norm of 
the difference between qj and quter— l) and the new 
density profile quer is set to this field qj. If the 
value selected for the angular velocity Sl^er is the 
corresponding to j = 0, then the value of (6£l)iter 
is set to half of (5il)( iter _i), otherwise it is left un- 
changed. 



Rest Mass: The Bernoulli constant C from Eq. i|34|) is 

adjusted to keep the desired value mo- 

Circulation: The spin parameter a is adjusted to keep 
the circulation at the desired value Coo- 

Tolerance Calculation: The new fields are compared 
with their respective values at the previous itera- 
tion step by means of the L2 norm. If the norm 
of the difference falls below some threshold value e 
for all the fields (i.e., index from 1 to 7), then the 
iteration is stopped. 



APPENDIX F: CALCULATION OF THE 
CIRCULATION USING PPN EXPANSIONS 

Using parametrized post-Newtonian expansions of the 
Einstein field equations we can estimate the relativistic 
circulation of fluid C around a closed path c inside the 
neutron star. We will consider the case of a binary in 
the inertial frame located at the system's center of mass. 
In this appendix, and (3 l will represent the fluid four- 
velocity and shift vector in the inertial frame. The binary 
is composed by a neutron star with parameters m s , V s , 
and r s , and a companion body with parameters m b , Vb, 
and r b . These parameters are the mass, orbital velocity, 
and coordinate position of the stellar center respectively. 
Additionally, we will assume that the velocity of the fluid 
behaves like a rigid rotator 

v = v s + [n s x (f-f s )Y , 

where Ct s is the spin angular velocity of the star. Note 
that we are ignoring any irrotational components of the 
velocity, like the term Va in Eq. I|29|l since, as it will be 
clear below, only the solenoidal part of the velocity will 
survive. 

To estimate the magnitude and scaling behavior of the 
fluid velocity we will take the components of the four- 
metric to be given by the series expansion 



.9oo 

90i 



-1 



2 4> N + 0(t 2 ) 
-4 Xl + 0(e 5 / 2 ) 



jij = (1 - 2 <M % + 0(e 2 ) . (Fl) 



Here e is the PPN expansion parameter (e ~ M/r ~ v 2 ) 
and the Newtonian potential (f>N and the potential \i are 
adapted from the expressions for the near-zone for point 
masses given by Will and Wiseman [H^ 



*>N = 4>s(r) - 



Xi 



-Ur)V l s 



n\ 

m b V b l 
\r- n\ 



(F2) 



where 4> s (r) is the self-gravitating Newtonian potential 
of the primary neutron star. We will ignore the effect of 
the spin of the bodies in the shift vector fl\ since it will 
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only contribute with terms that appear at higher order 
in e. 

The relativistic circulation is given by Eq. 1)30(1 



C(c) 



hujdx 1 



hv ■ dx 



V x(hv)-dS 



S(c 



[ft (V x v) + Vft x v\ ■ dS 



(F3) 



S(c) 



where we have used the fact that the line element along 
the closed path c is purely spatial der M = (0,dx l ), and 
Stokes theorem. The three-dimensional vector v has been 
defined by the spatial components Ui. We evaluate those 
components using Eq. I|C1|> and the metric ijFlfl 



= A + (1 - 3 <f> N )V l + 0(e 2 ) , 
where V — u l /u°. We then obtain 

v-P + V-Z <j) N V . (F4) 



With the definition of the PPN potentials l|F2J) . we obtain 
the following expressions for the curl of v (IF4I) 



4 m b r 
— — Wr-r 6 ) x V b 
\r - r b 3 



V x = 4V0 S x V s 

V x V = 20 s , 

and splitting the last term of the RHS of Eq. (|F4|) 

(f> N V = (j) N V s +(j} N [Q a x (f-f s )] 



(F5) 



= 

we get 
- 3V x (4> N V) = 



+ (/>jv[^ s x (r- r s )] 



- 3V0 s (r) x 
3 m;, 



f - n\ 



(r — n) xV s - 6(j) N fl s 



- W<t> N x [fi s x {f-f s )} . (F6) 

Let us now consider the particular case of an irrota- 
tional fluid. By definition, the vorticity of the fluid is 
null and thus the circulation of the fluid along any closed 
path c inside the star will be zero. Equation l|F3(l shows 
that this is only satisfied if 



h (V x v) + Vh x v = 



(F7) 



everywhere in the stellar volume, in particular the stel- 
lar center f = f s . The center is defined as the point 
of maximum density (or equivalently enthalpy), then 
V/i(f = f s ) = 0. Since ft. > 0, at the stellar center 
Eq. i|F7(l becomes V x v = 0. An expression for V x v 



can be constructed from Eqs. JF4J, (|F5|) . (|F6|) . Equating 
that expression to zero, setting f = f s , and solving for 
fi s gives 



^5 Ir^r's 



2\f s -fb\ 3 

(1-3 (M^))- 1 

m b 
2\f s -fb\ 3 



(f s - fb) x (3 V s - 4 Vb) 
(f s - f b ) x (3 V; - 4 H) (F8) 



where we use the fact that V</> s (r = f s ) — since 
4> s {f) oc (r — f s ) 2 in the neighborhood or f s . Note that 
the self-gravity terms in the metric cancel out of the final 
expression. Since we are working in the center of mass 
system and the objects are in circular orbit around the 
origin of coordinates, we have 

m b 

r, = d r, 

M 



Vs 

v b 



M 

Q x fs 
m s 

m b 



d f. 



V s , 



(F9) 



where d = \f s —f b \, f s = f s /r s , M = (m s +m b ), and f2 is 
the orbital angular velocity of the binary. Inserting the 
relations l|F9|l into Eq. (|F8|) gives 



n, = 



m b (2>m b + 4m s 

2 M d 
3m b + fi 



2 d 



(F10) 



where \i = (m b m s )/(m s + m b ) is the reduced mass. In 
the case of equal mass binaries like the ones studied on 
this paper, Eq. (|F10|) reduces to 

n s = • 

4 d 

Note that Eq. (|FT0)l is identical (to PPN order) to the 
formula for the precession of the spin of a star with 
respect to the inertial frame tied to distant galaxies. 
This corresponds to a combination of the Lense-Thirring 
(gravitomagnetic) and the geodetic effects |53| . 

The PPN lines shown on Fig. were calculated as 
{l s /Q — j m/d, where m is the gravitational mass of a 
star in isolation. 



APPENDIX G: TABLES 

This appendix summarizes the parameters correspond- 
ing to each circular orbit presented in this paper. Ta- 
bles III-X group the orbits belonging to a sequence of 
identical T = 2 binary polytropes with given compaction 
ratio (m/i?) oo and relativistic equatorial circulation C. 
The parameters tabulated are the coordinate separation 
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d, the total gravitational mass M, the binding energy 
BE = (M — Moo)/Mo, the total angular momentum 
J/Mq, the orbital angular frequency too, and the spin 
angular momentum Q s to . The total rest mass of the 
system with {m/R)^ = 0.14 (0.19) is M = 0.2938 
(0.3534). The critical value of equatorial circulation 
which corresponds to the case of single stars rotating at 
the mass-shedding limit is C crit = 1.96 (2.02) for the 
compaction ratio (m/R)^ — 0.14 (0.19). 
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TABLE III: (m/R)^ = 0.140, C/C cri t = -0.510. 



d 


M 


BE 




J/MS 


SI mo 




Sl s mo 




2.29 


0.2709 


-9.49 x 


10" 


-3 


0.7007 


1.87 x 


10" 


-2 


-1.44x 


10" 


-2 


2.41 


0.2710 


-8.98 x 


10~ 


-3 


0.7137 


1.74 x 


io- 


-2 


-1.48x 


io- 


-2 


2.57 


0.2713 


-8.24 x 


10~ 


-3 


0.7299 


1.58 x 


io- 


-2 


-1.51x 


io- 


-2 


3.06 


0.2716 


-6.97 x 


10~ 


-3 


0.7791 


1.28 x 


io- 


-2 


-1.56x 


io- 


-2 


3.43 


0.2719 


-6.04 x 


io- 


-3 


0.8187 


1.11 x 


io- 


-2 


-1.59x 


io- 


-2 


3.80 


0.2719 


-5.92 x 


io- 


-3 


0.8419 


9.26 x 


io- 


-3 


-1.53x 


io- 


-2 


4.46 


0.2723 


-4.79 x 


io- 


-3 


0.9061 


7.57 x 


io- 


-3 


-1.59x 


io- 


-2 


4.76 


0.2723 


-4.73 x 


io- 


-3 


0.9201 


6.78 x 


io- 


-3 


-1.54x 


io- 


-2 


4.97 


0.2724 


-4.31 x 


io- 


-3 


0.9430 


6.43 x 


io- 


-3 


-1.59x 


io- 


-2 


5.21 


0.2724 


-4.19 x 


10- 


-3 


0.9531 


5.95 x 


10- 


-3 


-1.60x 


10- 


-2 



TABLE IV: {m/R)oo = 0.140, C/C c nt = -0.255. 



d 


M 


BE 




J/M'Z 


too 




Sis too 




2.23 


0.2704 


-l.llx 


10" 


•i 


0.7273 


1.92x 


10" 


-'2 


-5.52x 


10" 


-a 


2.40 


0.2705 


-1.08x 


10" 


-2 


0.7475 


1.75x 


io- 


-2 


-5.72X 


10" 


-3 


2.57 


0.2707 


-l.Olx 


io- 


-2 


0.7688 


1.60x 


io- 


-2 


-5.90x 


10" 


-3 


3.26 


0.2712 


-8.26x 


io- 


-3 


0.8269 


1.16x 


io- 


-2 


-6.69 x 


10" 


-3 


3.94 


0.2716 


-6.97x 


io- 


-3 


0.8913 


8.92 x 


io- 


-3 


-7.09x 


10" 


-3 


4.63 


0.2719 


-6.03x 


io- 


-3 


0.9571 


7.04 x 


10" 


■3 


-7.10x 


10" 


-3 


5.14 


0.2720 


-5.57x 


10- 


-3 


0.9945 


6.08x 


10" 


-3 


-7.15x 


10" 


-3 



TABLE V: (m/R) x = 0.140, C/C cri t = 0.0. 



d 


M 


BE 




J/Mo 


fi mo 




Sis too 




1.89 


0.2698 


-1.31x 


10" 


-'2 


0.7494 


2.44x 


10" 


-'i 


3.29x 10" 


■3 


2.06 


0.2701 


-1.21x 


10" 


-2 


0.7689 


2.19x 


10" 


-2 


2.71 x 10 - 


■3 


2.38 


0.2705 


-1.07x 


10" 


-2 


0.8025 


1.82x 


10" 


-2 


1.93x 10" 


3 


2.74 


0.2708 


-9.93x 


10" 


-3 


0.8313 


1.48 x 


10" 


-2 


1.31x 10" 


■3 


3.49 


0.2713 


-8.15x 


10" 


-3 


0.8998 


1.06x 


10" 


-2 


0.73x 10" 


■3 


3.94 


0.2715 


-7.36x 


10" 


-3 


0.9421 


8.92x 


10" 


-3 


0.55x 10" 


-3 


4.71 


0.2719 


-6.06x 


10" 


-3 


1.0077 


6.98 x 


10" 


-3 


0.36x 10" 


-3 


5.46 


0.2721 


-5.50x 


10" 


-3 


1.0549 


5.56x 


10" 


-3 


0.24x 10" 


-3 


5.83 


0.2721 


-5.31 x 


10" 


-3 


1.0889 


5.11x 


10" 


-3 


0.22x 10" 


-3 



TABLE VI: {rn/R) x = 0.140, C/C cri t = 0.255. 



d 


M 


BE 




J/M£ 


Q r 






Sis TOO 




2.18 


0.2703 


-1.14x 


10" 


-'2 


0.8198 


2.02x 


10" 


•1 


8.12x 


10" 


-3 


2.52 


0.2706 


-1.03x 


10" 


-2 


0.8500 


1.65x 


10" 


-2 


7.87x 


10" 


-3 


2.74 


0.2709 


-9.56x 


10" 


-3 


0.8751 


1.49 x 


10" 


-2 


7.83 x 


10" 


-3 


3.26 


0.2712 


-8.29x 


10" 


-3 


0.9235 


1.18x 


10" 


-2 


7.51 x 


10" 


-3 


3.77 


0.2715 


-7.31 x 


10" 


-3 


0.9703 


9.69x 


10" 


-3 


7.38 x 


10" 


-3 


4.63 


0.2719 


-6.15x 


10" 


-3 


1.0413 


7.12 x 


10" 


-3 


7.27x 


10" 


-3 


5.14 


0.2720 


-5.63x 


10" 


-3 


1.0810 


6.08x 


10" 


-3 


7.25x 


10" 


-3 


5.66 


0.2721 


-5.24x 


10" 


-3 


1.1163 


5.19x 


10" 


•3 


7.19x 


10" 


-3 
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TABLE VII: (m/fl)oo = 0.140, C/C cri t = 0.510. 



d 


M 


BE 




J/M£ 


Q, mo 




£ls mo 




2.23 


0.2707 


i ni v 


1U 


-'1 


0.8716 


Z.KJo x 


1 A" 

1U 


-'2 


l.OOX 


1U 


-2 


2.44 


0.2710 


-9.40x 


10" 


-3 


0.8836 


1.77x 


10" 


-2 


1.33 x 


10" 


-2 


2.69 


0.2712 


-8.67x 


10" 


-3 


0.9072 


1.51x 


10" 


-2 


1.34 x 


10" 


-2 


2.77 


0.2712 


-8.41 x 


10" 


-3 


0.9159 


1.46 x 


10" 


-2 


1.35x 


10" 


-2 


3.03 


0.2715 


-7.63x 


10" 


-3 


0.9452 


1.27x 


10" 


-2 


1.36 x 


10" 


-2 


3.70 


0.2719 


-6.15x 


10" 


-3 


1.0128 


9.73 x 


10" 


-3 


1.38 x 


10" 


-2 


3.95 


0.2720 


-5.89x 


10" 


-3 


1.0275 


8.92x 


10" 


-3 


1.38 x 


10" 


-2 


4.79 


0.2723 


-4.89x 


10" 


-3 


1.0910 


6.77x 


10" 


-3 


1.41 x 


10" 


-2 


5.38 


0.2725 


-4.09 x 


10" 


-3 


1.1563 


5.71 x 


10" 


-3 


1.41 x 


10" 


-2 


5.97 


0.2726 


-3.90x 


10" 


-3 


1.1779 


4.81 x 


10" 


-3 


1.42 x 


10" 


-2 



TABLE VIII: {m/R) x = 0.190, C/C cri t = -0.124. 



d 


M 


BE 




J/MS 


f2 mo 




Q s mo 




2.69 


0.3186 


-l.lOx 


10" 


•i 


0.7331 


1.93x 


10" 


-'2 


-4.08x 10" 


-3 


3.10 


0.3190 


-9.92x 


10" 


-3 


0.7608 


1.58x 


10" 


-2 


-4.66x 10" 


-3 


3.28 


0.3192 


-9.34x 


10" 


•3 


0.7838 


1.48 x 


10" 


-2 


-4.78x 10" 


-3 


3.51 


0.3193 


-8.96x 


10" 


-3 


0.7956 


1.34x 


10" 


-2 


-4.92x 10" 


-3 


4.03 


0.3197 


-7.97x 


10" 


-3 


0.8358 


l.llx 


10" 


-2 


-5.24x 10" 


-3 


4.29 


0.3198 


-7.59x 


10" 


-3 


0.8516 


l.Olx 


10" 


-2 


-5.39x 10" 


-3 


5.04 


0.3202 


-6.61x 


10" 


-3 


0.9089 


8.07x 


10" 


-3 


-5.50x 10" 


-3 



TABLE IX: (m/R) x = 0.190, C/C crlt = 0.0. 



d 


M 


BE 




J/Mo 


Q mo 




O s mo 




2.52 


0.3182 


-1.19x 


10" 


-2 


0.7415 


2.11x 


10" 


-2 


2.34x 


10" 


-3 


2.86 


0.3187 


-1.08x 


10" 


-2 


0.7615 


1.77x 


10" 


-2 


1.69 x 


10" 


-3 


3.03 


0.3189 


-1.02x 


10" 


-2 


0.7799 


1.65x 


10" 


-2 


1.52x 


10" 


-3 


3.45 


0.3192 


-9.34x 


10" 


-3 


0.8036 


1.36x 


10" 


-2 


1.07x 


10" 


-3 


3.70 


0.3194 


-8.69 x 


10" 


-3 


0.8308 


1.25x 


10" 


-2 


0.94x 


10" 


-3 


4.37 


0.3198 


-7.63x 


10" 


-3 


0.8772 


9.86x 


10" 


-3 


0.62x 


10" 


-3 


4.96 


0.3201 


-6.87x 


10" 


-3 


0.9193 


8.26x 


10" 


-3 


0.46 x 


10" 


-3 


5.72 


0.3202 


-6.42 x 


10" 


-3 


0.9481 


6.58x 


10" 


•3 


0.32x 


10" 


-3 



TABLE X: (m/R)^ = 0.190, C/C cri t = 0.247. 



d 


M 


BE 




J/Mo 


Q mo 




O s mo 




2.64 


0.3185 


-l.llx 


10" 


-'2 


0.7879 


2.00 x 


10" 


-2 


1.32x 


10" 


-2 


2.77 


0.3186 


-1.07x 


10" 


-2 


0.7977 


1.85x 


10" 


-2 


1.31x 


10" 


-2 


3.12 


0.3189 


-9.79x 


10" 


-3 


0.8224 


1.55x 


10" 


-2 


1.26x 


10" 


-2 


3.57 


0.3193 


-8.67x 


10" 


-3 


0.8642 


1.29x 


10" 


-2 


1.29x 


10" 


-2 


4.09 


0.3195 


-8.07x 


10" 


-3 


0.8899 


1.09x 


10" 


-2 


1.30 x 


10" 


-2 


4.36 


0.3197 


-7.49 x 


10" 


-3 


0.9189 


9.61x 


10" 


-3 


1.30 x 


10" 


-2 


5.80 


0.3204 


-5.86x 


10" 


-3 


0.9954 


6.48 x 


10" 
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